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Abstract 

We  consider  the  problem,  of  determining  functions  of  an  image  of  an  object  that  are  insensitive 
to  illumination  changes.  We  first  show  that  for  an  object  with  Lambertian  reflectance  there  are 
no  discriminative  functions  that  are  invariant  to  illumination.  We  do  this  by  showing  that  given 
any  two  images,  one  can  construct  a  single  Lambertian  object  that  can  produce  both  images  under 
two  very  simple  lighting  conditions.  This  result  leads  us  to  adopt  a  probabilistic  approach  in  which 
we  analytically  determine  a  probability  distribution  for  the  image  gradient  as  a  function  of  the 
surface ’s  geometry  and  reflectance.  Our  distribution  reveals  that  the  direction  of  the  image  gradient 
is  insensitive  to  changes  in  illumination  direction.  We  verify  this  empirically  by  constructing  a 
distribution  for  the  image  gradient  from  more  than  20  million  samples  of  gradients  in  a  database  of 
1,280  images  of  20  inanimate  objects  taken  under  varying  lighting  conditions.  Using  this  distribution, 
we  develop  an  illumination  insensitive  measure  of  image  comparison  and  test  it  on  the  problem  of 
face  recognition. 


‘This  paper  is  based  on  “Comparing  Images  Under  Variable  Illumination,”  by  Jacobs,  Belhumeur,  and  Basri,  which 
appeared  in  the  IEEE  Conference  on  Computer  Vision  and  Pattern  Recognition,  June  1998,  and  on  “In  Search  of 
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1  Introduction 


Changes  in  viewpoint  and  illumination  can  dramatically  alter  the  appearance  of  an  object.  To 
generalize  effectively,  image-based  recognition  systems  must  use  methods  of  image  comparison  that 
work  in  spite  of  these  changes.  We  focus  here  solely  on  changes  in  illumination  and  ask:  Are  there 
discriminative  illumination  invariants?  If  not,  are  there  local  image  measurements  that  are  at  least 
insensitive  to  illumination  changes? 

We  show  that  even  for  objects  with  Lambertian  reflectance  [28],  there  are  no  discriminative 
functions  of  images  of  objects  that  are  invariant  to  illumination.  This  differs  from  earlier  findings  in 
that  we  do  not  assume  a  homogeneous  BRDF  [27,  5],  coplanarity  [36],  or  consider  invariants  based 
on  multiple  images  [48].  To  do  this,  we  show  that  for  any  two  images  —  whether  or  not  they  are 
of  the  same  object  —  there  is  always  a  family  of  surfaces,  albedo  patterns,  and  light  sources  that 
could  have  produced  them. 

This  result  suggests  that  in  comparing  images  for  recognition,  alignment,  or  tracking,  one  must 
resort  to  probabilistic  measures  of  comparison.  (Here  we  do  not  assume  we  have  access  to  multiple 
training  images  of  the  same  object  under  varying  illumination.  If  we  did,  we  could  extract  invariants 
[48]  or  construct  representations  for  modeling  the  illumination  variation  [43,  20,  3,  18].)  We  develop 
a  measure  of  image  comparison  by  considering  illumination  to  be  a  random  variable  which  gives  rise 
to  the  apparent  randomness  in  local  image  measurements.  Specifically,  we  derive  the  probability 
distribution  of  the  image  gradient  of  a  point  on  a  surface  as  determined  by  the  differential  geometric 
and  reflectance  properties  at  that  point.  Our  distribution  reveals  that  the  direction  of  the  image 
gradient  (which  is  perpendicular  to  the  flow  field  [5])  is  insensitive  to  changes  in  illumination  direc¬ 
tion.  Using  the  image  gradient  distribution,  we  then  develop  an  illumination  insensitive  measure  of 
image  comparison. 

Finally,  we  verify  our  qualitative  theoretical  arguments  by  empirically  constructing  both  our 
distribution  for  the  image  gradient  and  our  illumination  insensitive  measure  of  image  comparison. 
We  do  this  using  a  database  of  1,280  images  of  20  objects  taken  under  varying  illumination  direction. 
We  then  test  our  illumination  insensitive  measure  on  the  problem  of  face  recognition,  and  compare 
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our  performance  on  450  images  of  10  individuals  to  that  of  other  existing  methods. 

2  Background 

In  recent  years  there  has  been  much  work  on  object  recognition  by  image  comparison.  In  these 
methods,  an  object  is  not  described  in  terms  of  its  3-D  properties,  but  rather  in  terms  of  the  2-D 
images  that  it  produces.  One  approach  to  appearance-based  recognition  is  to  sample  an  object’s 
possible  images,  and  then  to  compare,  in  a  lower  dimensional  image  subspace,  a  novel  image  to  the 
set  of  sampled  images,  using  pattern  recognition  techniques  such  as  nearest  neighbors.  Turk  and 
Pentland  [44]  (inspired  by  the  work  of  Kirby  and  Sirovich  [26])  suggest  such  an  approach,  using 
principal  component  analysis  to  compactly  represent  the  training  images.  Murase  and  Nayar  [35] 
suggest  accounting  for  lighting  and  viewpoint  variation  with  such  an  approach  by  explicitly  sampling 
images  of  the  object  under  all  possible  viewing  conditions.  The  advantage  of  these  methods  is  that 
they  do  not  have  to  derive  3-D  structure,  or  to  extrapolate  to  all  possible  images  of  an  object  based 
on  a  small  set  of  samples.  The  disadvantage  is  that  the  set  of  images  an  object  can  produce  is 
extremely  large,  since  there  are  so  many  ways  that  viewing  conditions  can  vary.  Murase  and  Nayar, 
for  example,  are  only  able  to  deal  with  one  rotational  degree  of  freedom  in  the  object,  and  one 
positional  degree  of  freedom  in  a  point  light  source.  In  principal  their  approach  can  handle  greater 
variability,  but  the  number  of  images  an  object  produces  grows  exponentially  with  the  number  of 
degrees  of  freedom  in  the  viewing  conditions. 

To  build  flexible  recognition  systems  that  account  for  all  the  variability  of  the  image  formation 
process,  it  seems  necessary  to  find  ways  to  generalize  from  images  that  are  few  relative  to  this 
variability.  The  most  powerful  form  of  generalization  is  based  on  invariants.  An  invariant  is  a 
property  of  an  object  that  shows  up  in  every  image  of  that  object.  Geometric  invariants  have  been 
used  for  the  recognition  of  restricted  classes  of  objects,  (e.g.,  planar  [29]  or  symmetric  [16])  that  have 
identifiable  local  features.  But  it  has  been  shown  that  invariants  do  not  exist  for  general  classes  of  3- 
D  objects  ([7,  9,  33]).  Photometric  invariants  have  also  been  used  for  planar  objects  ([36,  42]).  Moses 
and  Ullman  have  shown  that  photometric  invariants  do  not  exist,  even  for  Lambertian  objects,  if  the 
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“objects”  are  allowed  to  consist  of  mosaics  of  possibly  disconnected  planar  patches.  In  this  paper, 
we  show  the  stronger  result  that  discriminative  invariants  do  not  exist  even  when  these  surface 
normals  are  required  to  form  surfaces.  Both  geometric  and  photometric  invariants  can  also  exist 
when  one  has  access  to  multiple  images  of  an  object  at  recognition  time.  Such  techniques  are  used 
for  structure-from-motion  (e.g.,  [13])  or  photometric  stereo  (e.g.,  [48,  10,  37]). 

In  the  absence  of  invariants  for  recognition,  one  can  attempt  to  use  a  small  number  of  images 
of  an  object  to  characterize  the  complete  set  of  images  that  the  object  might  produce.  Ullman 
and  Basri  [45]  take  this  approach,  showing  that  in  many  cases  each  of  the  2-D  images  produced 
by  an  object  from  different  viewpoints  is  a  linear  combination  of  a  small  number  of  basis  images. 
The  spirit  of  this  work  is  to  use  images  of  an  object  to  predict  what  other  images  the  object  might 
produce,  without  ever  explicitly  reconstructing  the  object’s  3-D  properties.  Jacobs  [23]  determines 
the  most  compact  possible  representations  of  these  sets  of  images,  while  Faugeras  and  Robert  [14] 
presents  methods  for  extrapolating  the  images  that  an  object  produces  under  perspective  projection. 
Shashua  [43]  and  Moses  [31]  take  an  approach  similar  to  that  of  Ullman  and  Basri  to  show  that  the 
set  of  intensity  images  produced  by  an  object  under  variable  lighting  conditions  can  be  predicted 
linearly  from  a  small  number  of  basis  images.  Belhumeur  and  Kriegman  [3]  show  how  to  extend 
this  work  to  account  for  multiple  or  diffuse  light  sources,  and  self-shadowing.  In  all  this  work,  it 
is  assumed  that  training  images  provide  sufficient  prior  information  to  completely  characterize  the 
entire  set  of  images  that  an  object  can  produce.  This  information  may  be  strictly  weaker  than 
having  a  precise  3-D  description  of  the  object,  because  various  approaches  may  not  attempt  to 
recover  elements  of  the  structures  that  are  altered  by  a  3-D  affine  transformation  (e.g.,  [45,  43])  or  a 
bas-relief  transformation  (e.g.,  [4]).  However,  it  is  clear  that  much  of  the  structure  of  an  object  must 
be  recovered  in  order  to  predict  all  images  that  it  can  produce,  and  that  this  may  be  difficult  to  do 
based  on  a  small  set  of  unconstrained  training  images.  For  example,  Jacobs,  Belhumeur,  and  Basri 
[24]  show  that  under  arbitrary,  diffuse  lighting  conditions,  it  is  not  possible  to  exactly  reconstruct 
the  albedos  of  an  object,  even  when  its  structure  is  precisely  known. 

Finally,  one  may  wish  to  perform  recognition  under  variable  viewing  conditions  when  there 
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is  simply  not  enough  prior  information  to  fully  characterize  the  effect  of  this  variation.  In  this 
case,  one  may  wish  to  use  image  comparision  methods  that  are  insensitive  to  changes  in  viewing 
condition,  even  if  they  are  not  completely  invariant  to  changes.  These  are  sometimes  called  quasi¬ 
invariants.  A  number  of  geometric  quasi-invariants  have  been  proposed  (e.g.,  [7,  51]);  they  are 
sometimes  relaxed  versions  of  invariant  properties.  Related  photometric  methods  have  long  been 
used;  one  justification  for  edge-based  recognition  methods  is  that  for  polyhedral  objects  edge  position 
is  insensitive  to  illumination  and  viewpoint  changes.  (However  for  smooth  surfaces,  edge  position 
is  sensitive  to  illumination  changes.)  More  recently,  Rao  and  Ballard  [39],  for  example,  recognize 
objects  using  descriptors  built  from  the  multi-scale  output  of  filters.  These  are  insensitive  to  modest 
changes  in  viewpoint  and  illumination  (for  insensitive  comparisons  of  color  see,  e.g.,  [17]).  Irani  and 
Anandan  [22]  suggest  comparing  images  using  the  squared  output  of  directional  derivative  filters  to 
gain  insensitivity  to  changes  in  sensing  modality. 

After  proving  that  illumination  invariant  properties  do  not  exist  for  3-D  objects,  we  also  develop 
a  probabilistic  semi- invariant  measure,  namely,  the  direction  of  the  image  gradient,  that  is  insensitive 
to  lighting  variation.  Our  method  is  based  on  a  geometric  analysis  of  the  effects  of  lighting  variation 
on  the  image  gradients  produced  by  a  3-D  object.  We  compare  this  to  empirical  data,  and  produce 
and  test  a  new  image  comparison  method. 

3  Do  Illumination  Invariants  Exist? 

Given  two  images,  are  they  created  by  two  distinct  objects,  or  the  same  object  under  different 
illuminations?  While  it  seems  possible  that  one  can  always,  or  at  least  sometimes,  distinguish  with 
certainty  between  these  two  scenarios,  in  this  section  we  show  the  contrary  is  true. 

To  analyze  this  question,  we  study  the  existence  of  discriminative  illumination  invariants:  func¬ 
tions  of  an  image  which  are  invariant  to  illumination  on  the  object  but  vary  with  the  object  identity. 
Formally,  let  O  be  some  set  of  rigid  objects  including  their  optical  properties;  S  be  certain  illumi¬ 
nation  conditions;  and  I  be  the  set  of  all  images.  Function  Q  :  O  x  S  — >  I  gives  the  image  I  E  T  of 
object  ogO  under  illumination  s  E  5,  i.e. ,  I  =  Q(o,  s). 
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We  adopt  the  following  definitions: 


Definition  1  A  function  //,  on  I  is  invariant  to  illumination  •<=>  p(Q(o,s))  =  / x(Q(o,l )),  Vs,  l  E 
5,o  E  O. 

Definition  2  An  illumination  invariant  n  is  nondiscriminative  for  object  set  O  •<=>  n(I)  =  n(J), 
VI  J,  where  I  and  J  are  in  the  range  of  Q.  An  illumination  invariant  is  discriminative  iff  it  is 
not  nondiscriminative. 

This  definition  implies  /i  does  not  depend  on  o  for  nondiscriminative  invariants. 

Lemma  3.1  There  are  no  discriminative  illumination  invariants  for  O  if  for  any  two  images  I  and 
J  in  the  range  of  Q,  there  is  always  an  object  o  E  O  which,  under  some  pair  of  lighting  conditions 
in  S,  could  have  produced  both  image  I  and  J. 

Proof.  The  proof  follows  immediately  from  the  above  definitions.  □ 

The  existence  of  discriminative  invariants  depends  on  the  specific  set  of  illumination  and  objects. 
It  is  obvious  that  for  a  large  enough  class  of  illuminations,  there  are  no  such  invariants.  For  example, 
a  movie  projector  can  create  arbitrary  images  by  projecting  carefully  designed  patterns  on  any  given 
object  with  nonzero  reflectance.  The  larger  the  set  of  illumination  conditions  we  admit,  the  smaller 
the  set  of  invariants  will  be.  Likewise,  for  purely  specular  convex  surfaces,  (e.g.,  mirrors)  it  is  again 
obvious  that  there  are  no  discriminative  illumination  invariants.  However  these  cases  are  clearly 
extremes  in  illumination  and  reflectance.  Next,  we  consider  in  detail  the  case  of  Lambertian  surfaces 
illuminated  by  point  light  sources  distant  from  the  object.  What  is  surprising  is  that  even  for  this 
simple  case  there  are  no  discriminative  illumination  invariants.  We  first  prove  the  claim  ignoring  the 
effects  of  interreflection,  then  show  that  our  result  holds  when  interreflection  is  taken  into  account. 
Moreover,  our  simulations  show  that  their  additive  effect  on  the  image  intensity  is  often  negligible 
when  considered  in  the  context  of  this  problem  [8].  The  fact  that,  statistically  the  typical  albedo  of 


6 


an  object  is  less  than  0.15  [40,  41],  lends  empirical  support  for  the  technique,  letting  the  albedo  of 
the  object  approach  zero,  employed  in  the  proof  of  Theorem  3.2  below. 

Theorem  3.1  Discounting  interreflection,  given  two  arbitrary  image  functions  I  and  J,  and  two 
arbitrary  point  light  sources  s  and  l  at  infinity,  and  if  the  projection  on  the  mage  plane  of  the  sum 
s  +  l  is  nonzero,  there  exists  a  family  of  smooth  Lambertian  surface  f  such  that  I  is  the  image  of  f 
under  s  and  J  is  that  of  f  under  l. 

Proof.  This  is  just  a  recapitulation  of  the  Lemmas  3.5  and  3.2  that  are  to  follow.  □ 

Corollary  3.1  Discounting  interreflection,  all  illumination  invariants  for  objects  with  Lambertian 
reflectance  under  point  light  sources  at  infinity  are  nondiscrim, inative. 

Proof.  Let  ft  be  an  illumination  invariant.  Given  two  arbitrary  images  7,  J  in  the  range  of  Q, 
by  Theorem  3.1,  there  always  exists  a  Lambertian  surface  of  an  object  o  and  two  light  sources  at 
infinity  s  and  l,  such  that  I  and  J  are  images  of  o  under  s  and  l,  respectively,  or  I  =  Q(o,s)  and 
J  =  Q(o,l).  By  Lemma  3.1,  //,  must  be  a  nondiscriminative  invariant.  □ 

Note  that  for  the  above  propositions  we  require  objects  to  be  composed  of  surfaces,  not  freely 
floating  planar  facets  in  space  as  in  [32], 

We  prove  the  much  stronger  Theorem  3.1  rather  than  the  weaker  Corollary  3.1  directly,  since 
some  researchers  believe  humans  are  able  to  determine  the  direction  of  the  light  source  when  viewing 
an  image  [38,  6,  21,  25,  30,  46,  49,  50,  52],  We  show  even  under  this  more  stringent  condition,  there 
is  still  no  way  to  tell  with  certainty  if  the  given  two  images  are  generated  from  two  different  surfaces. 

If  we  assume  in  addition  that  the  invariant  function  is  continuous  with  respect  to  the  image 
functions,  the  above  result  can  be  extended  to  Lambertian  surfaces  with  interreflection. 

Theorem  3.2  Taking  interreflection  into  account,  all  continuous  illumination  invariants  for  objects 
with  Lambertian  reflectance  under  point  light  sources  at  infinity  are  nondiscriminative. 
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Proof.  See  appendix. 


□ 


In  this  section  we  outline  the  proof  of  Theorem  3.1.  The  basic  strategy  of  our  proof  is  to 
write  two  equations  that  describe  the  formation  of  the  images,  with  the  object’s  shape  and  albedo 
as  unknowns.  After  eliminating  albedo,  we  arrive  at  a  single  first  order  linear  partial  differential 
equation  in  shape.  We  then  show  that  this  PDE  always  has  a  unique  solution,  up  to  a  set  of  initial 
conditions  that  determine  a  family  of  solutions.  Using  this  object  shape,  albedo  is  uniquely  solved 
for  linearly,  up  to  a  scale  factor.  Finally,  we  must  show  that  the  solutions  we  arrive  at  do  not  contain 
shadows,  which  justifies  our  use  of  a  simple  linear  description  of  image  formation. 

We  begin  the  outline  of  our  proof  of  Theorem  3.1  with  the  following  setup.  Set  the  camera 
optical  axis  as  the  z-axis.  An  object  o  is  viewed  from  the  direction  (0,0,1).  The  set  O  (more 
accurately  the  visible  surface  of  O)  is  defined  to  be  the  set  of  pairs  (/,  a),  where  /  is  the  graph 
( x,y,z  =  f(x,y))  of  a  smooth  function  z  =  f(x,y ),  and  a(x,y)  is  a  nonnegative  function  called 
albedo,  both  defined  on  the  prescribed  unit  square.  Here,  the  set  of  point  light  sources  at  infinity 
S  is  represented  by  the  set  of  3-D  vectors  s  =  (sxi sy.  sz)  in  the  opposite  direction  of  the  light  rays, 
with  P|  equal  to  the  powers  of  the  sources.  There  are  two  regions  on  the  surface  (x,  y ,  f{x,  y))\  one 
that  is  illuminated  by  the  light  source,  and  the  other  called  shadow  region  (SR  hereafter)  that  is 
either  facing  away  and  creates  so  called  attached  shadows,  or  obscured  from  the  light  source  by  other 
parts  of  the  surface.  Suppose  an  image  I  is  generated  by  /  under  point  light  sources  s  =  (sx,  sy,  sz). 
Excluding  interreflection,  the  image  intensity  I  is  described  by 


f  a{x,y)  s-h{x,y)  if  (x,  y,  f(x,  y))  £  SR 
\  0  otherwise 


(1) 


where  a(x,y)  is  the  albedo  of  object  o,  n  =  (  — 1|,  —  1)  is  a  normal  vector  to  the  surface  /  of  o, 

and  n  is  the  unit  vector  for  n.  Henceforth,  we  will  adopt  the  following  convention:  if  r  is  a  vector, 
then  r  =  r/||r  ||. 

There  is,  in  principle,  no  more  restriction  other  than  nonnegativity  for  a  function  to  describe 
the  image  of  some  object.  By  painting  the  albedo  proportional  to  the  given  image  function  on  a 
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flat  surface,  and  placing  the  surface  under  any  point  light  source  with  appropriate  power,  one  can 
always  obtain  any  prescribed  image.  However,  for  convenience  of  later  discussion,  we  define 

Definition  3  I  is  called  an  image  function,  or  I  E  T,  if  I  is  C1  and  positive  on  a  compact  set.  C1 
on  the  boundary  means  there  is  a  tangent  plane  approached  by  interior  and  boundary  points. 

The  restriction  is  not  severe  in  that  any  real  valued  piecewise  continuous  function  on  a  compact 
set  can  be  approached  pointwise  by  a  sequence  of  continuous  functions  which  can  in  turn  be  ap¬ 
proached  uniformly  by  a  sequence  of  analytic  functions  (Stone- Weierstrass  approximation  theorem). 
We  may  directly  introduce  discontinuities.  However,  its  existence  and  the  specific  geometry  of  the 
set  on  which  /  is  discontinuous  will  greatly  complicate  the  ensuing  technical  investigation  without 
altering  the  essential  result1. 

Our  goal  is  to  show  that  for  any  two  I,J  E  X,  there  always  exists  a  smooth  solution  on  the  unit 
square  to  the  partial  differential  equations  (PDE’s)  given  in  Eq.’s  2  and  3,  which  in  turn  leads  to 
Theorem  3.1. 

We  construct  a  surface  according  to  the  local  differential  description,  or  the  upper  part  of  Eq.  1, 
as  if  there  is  no  effect  of  cast  shadows. 


l{x,y) 

=  a{x,y)  s-n{x,y), 

(2) 

J  (x)  y) 

=  a(x,  y)  l  ■  n(x,  y) 

(3) 

We  then  show  that  if  the  boundary  conditions  are  chosen  appropriately,  the  surface  thus  con¬ 
structed  actually  casts  no  shadows,  and  therefore  renders  the  lower  condition  of  Eq.  1  moot: 

Lemma  3.2  Let  H  be  a  compact  subset  o/]R2.  Suppose  I,  J  E  Z(fi)  (image  functions  on  Q).  There 
is  a  family  of  C 1  surfaces  f  that  satisfies  Eq.  2  and  3  has  no  attached  or  cast  shadows. 

To  prove  the  above  lemma,  we  make  the  following  construction  which  leads  to  three  additional 

1If  there  are  discontinuities,  if  the  set  of  discontinuities  obey  certain  smoothness  conditions,  the  set  of  (/,  s)  that 
admits  continuous  surfaces  may  not  be  the  whole  6-D  space,  but  will  be  of  nonzero  measure 
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lemmas.  (The  trusting  reader  may  want  to  skip  ahead  to  Section  4.)  All  the  proofs  are  relegated  to 
the  Appendix. 


We  cross  multiply  the  two  sides  of  Eq.  2  with  those  of  Eq.  3.  We  then  divide  the  resulting 
equation  by  a  (a  cannot  be  zero  as  neither  I  nor  J  assumes  zero  value  anywhere)  to  obtain  the 
following  first  order  linear  partial  differential  equation: 

(if-  Js)  ■  n  =  0.  (4) 


(Wolff  and  Angelopoulou  [47]  and  Fan  and  Wolff  [12]  also  derive  this  equation,  and  then  use  it  for 
stereo  matching  and  photometric  stereo).  Once  this  equation  is  solved,  we  can  substitute  the  result 
into  our  original  equation  to  obtain  a  linear  equation  in  albedo.  When  the  magnitude  of  the  light 
is  unknown,  this  equation  determines  albedo  uniquely,  up  to  a  scale  factor.  To  obtain  valid  albedos 
we  must  scale  the  magnitude  of  the  light  so  that  all  albedos  are  less  than  one. 

To  solve  Equation  4  we  use  the  method  of  characteristic  curves.  In  this  method,  one  performs 
a  change  of  variables  to  obtain  a  PDE  in  one  variable.  Solving  this  PDE  tells  us  the  surface  height 
along  a  curve  on  the  object.  The  complete  equation  will  have  a  unique  solution,  up  to  an  initial 
condition,  if  every  point  is  on  a  single  characteristic  curves,  so  that  solving  these  curves  separately 
provides  exactly  one  height  for  every  point  on  the  surface. 

As  a  simple  example  of  this  method,  consider  the  case  of  l  =  (0,  0, 1)  and  s  =  (1,  0,  0).  Then  we 
have 


ox 

Our  equation  reduces  to  one  with  a  partial  derivative.  The  characteristic  curves  in  this  case  are 
horizontal  lines  across  the  image.  The  value  of  f(x,y)  is  given  by 


f(x,y)  =  f{0,y)+l  V |  dw 

Jo  J(w,y) 

where  we  have  no  obvious  source  of  knowledge  available  to  provide  the  initial  condition  /(0,y). 

The  above  construction  provides  one  pair  of  point  light  sources  that  when  provided  with  an 
appropriate  boundary  condition  gives  the  required  smooth  surface.  However,  here  we  intend  to 
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show  the  existence  of  a  surface  even  when  the  point  light  sources  are  given  for  each  of  the  two 
images. 

Generally,  the  characteristic  curves  r(t)  of  Eq.  4  satisfy 

dr  -» 

—  =Il-  Js.  (5) 

d  t  v  ' 

Obviously,  we  only  need  to  solve  for  the  characteristic  curves  in  the  XY  plane.  The  z-component 
can  be  solved  by  simply  integrating  the  z-component  of  the  right  hand  side  above  along  the  XY 
trajectory.  The  surface  is  constructed  by  weaving  the  characteristic  curves  together  by  choosing 
continuous  boundary  values  on  part  of  the  boundary  of  the  unit  square. 

We  claim  that  for  a  vector  field  uniformly  bounded  from  below,  as  is  the  case  here,  there  is  a 
characteristic  curve  in  the  XY-plane  through  each  point  in  D 2  with  both  ends  lying  on  the  boundary 
dD2. 


Lemma  3.3  Let  f l  be  a  compact  subset  of  JR2.  Suppose  I,  J  E  T(Ll)  (image  functions  on  Ll),  and 
the  projection  on  the  image  plane  of  the  sum  s  +  l  is  nonzero.  Let  p  denote  points  in  the  XY-plane. 
Through  an  arbitrary  point  q  in  f l,  there  exists  a  unique  (q  E  R  and  a  unique  characteristic  curve 
projection  on  the  XY-plane  p(t),  t  E  [0,  £5]  in  ^  homeomorphic  to  [0,  (q],  such  that  p(0),p((q)  E  dLl. 

Lemma  3.3  shows  the  global  existence  of  characteristic  curves  through  any  point  in  the  unit  square. 

Lemma  3.4  Let  the  hypothesis  of  Lemma  3.3  be  satisfied.  The  XY-plane  characteristic  curve  is  C1 
with  respect  to  t  and  its  initial  point. 

The  uniqueness  of  the  characteristic  curves  provided  by  Lemma  3.4  removes  the  ambiguity,  enabling 
the  later  construction  of  a  smooth  surface  instead  of  a  multivalued  function.  It  also  guarantees  the 
needed  smoothness  for  the  characteristic  curves,  although  not  yet  for  the  surface. 

The  smooth  surface  can  then  be  fabricated  by  carefully  assigning  smooth  boundary  conditions. 

Lemma  3.5  Let  the  hypothesis  of  Lemma  3.3  be  satisfied.  There  is  a  C 1  surface  f  on  f l  satisfying 
Eq.  4. 
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The  above  construction  from  the  local  description  would  not  satisfy  Eq.  1  if  the  surfaces  cast 
shadows  onto  themselves.  Our  task  would  be  made  much  simpler  if  there  is  a  surface  without 
attached  or  cast  shadows.  Fortunately,  Lemma  3.2  assures  that  there  indeed  exist  such  surfaces 
provided  that  the  boundary  condition  is  judiciously  chosen. 

To  summarize,  we  have  shown  that  given  any  two  images,  and  any  two  point  light  sources,  we 
can  construct  a  Lambertian  surface  that  will  generate  each  image,  in  the  absence  of  interreflections. 

In  the  appendix,  we  show  that  the  effects  of  interreflections  can  be  made  arbitrarily  small.  This,  in 
turn,  shows  that  no  discriminative  illumination  invariants  exist. 

4  A  Probabilistic  Model 

In  the  previous  section,  we  concluded  that  for  Lambertian  and  purely  specular  surfaces  there  are 
no  illumination  invariants.  Thus,  we  must  settle  for  something  less:  a  weaker,  probabilistic  form 
of  invariance.  To  this  end,  we  will  show  that  even  if  the  direction  of  the  light  source  is  uniformly 
distributed,  the  direction  of  the  image  gradient  is  not.  We  first  analytically  construct  the  probability 
distribution  for  the  image  gradient  as  a  function  of  the  surface  curvature  and  reflectance.  We  will 
verify  this  empirically  by  constructing  a  distribution  for  the  image  gradient  from  a  database  of  real 
images. 

Suppose  (x,y,  f  (x,y))  is  a  smooth,  i.e.  Ck  or  analytic,  surface.  Let  (u,v)  be  coordinates  on  the 
surface  such  that  (x(u,  •),  y(u,  •),  f[x(u,  ■ ),y{u ,  •)])  as  a  function  of  u  only  and  (x(-,  v),y(-,  v),  f[x(-,  v),y(-,v)]) 
as  a  function  of  v  only  are  lines  of  curvature,  with  u,  v  being  the  length  of  the  respective  lines.  For 
subsequent  development,  let  k.u  and  kv  be  the  principal  curvatures  of  the  surface  in  principal  direc¬ 
tions  u  and  v,  respectively.  We  set  up  a  local  Cartesian  coordinate  system  by  adopting  the  tangents 
of  the  lines  of  curvature  along  with  the  surface  normal  as  the  third  axis  and  call  this  the  u-v-n 
coordinate  system. 

The  BRDF  (bidirectional  reflectance  distribution  function)  a  is  a  function  of  (u,  v),  the  direction 
of  the  incident  light  i.  and  the  direction  of  outgoing  light  6  in  this  local  coordinate  system.  The 
radiance  from  a  point  on  a  smooth  surface  in  the  camera  direction  c  (in  the  u-v-n  coordinate  system) 
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is  then 


L(u,v,c)  =  /  a[u,  v,  i(u,  v,  s),  o(u,  «,£)]  h(u,  v)  ■  s  ds 

Jn 

=  h  ■  as  dJ  (6) 

J  n 

where  D  denotes  the  solid  angle  of  light  seen  at  the  point  of  concern. 

We  analyze  the  influence  of  the  differential  geometric  and  reflectance  properties  by  examining 
the  scene  radiance  under  a  single  light  source  at  infinity.  The  scene  radiance  in  Eq.  6  becomes 

L(u,  v,  c)  =  a[u ,  v,i(u ,  v,  s),  d(u,  v,  c)]  s-  n(u,v).  (7) 

Assume  the  light  sources  are  spherically  symmetrically  distributed.  Suppose  we  have  a  patch  of 
a  Lambertian  surface  with  constant  albedo  and  principal  curvatures  \ku\  /  \kv\.  Then  the  gradient 
of  the  scene  radiance  most  likely  lies  in  the  direction  in  which  the  magnitude  of  the  curvature 
is  maximal.  Note  that  in  the  special  case  of  cylindrical  objects,  e.g.,  coffee  mugs,  animal  limbs, 
telephone  poles,  etc.,  the  image  intensity  gradient  lies  orthogonal  to  the  cylinder  axis  and  tangent 
to  the  surface,  and  the  isophotes  run  parallel  to  the  cylinder  axis.  As  a  second  example,  let  us 
consider  a  planar  patch  of  surface  with  nonhomogeneous  BRDF.  For  any  BRDF,  the  direction  of 
the  gradient  of  the  scene  radiance  always  points  in  the  direction  of  the  spatial  gradient  of  the 
BRDF.  These  observations  suggest  that  if  the  light  source  directions  are  distributed  uniformly, 
the  distribution  of  the  gradient  of  scene  radiance  will  be  biased  by  the  underlying  geometry  and 
reflectance. 

We  shall  expound  these  observations  in  two  steps.  First,  we  derive  the  relation  between  the 
gradient  of  the  scene  radiance  VL  and  the  surface’s  local  geometry  and  reflectance.  Second,  we 
impose  a  spherically  symmetric  probability  distribution  on  the  light  source  s  and  determine  the 
resulting  distribution  on  the  gradient  of  the  scene  radiance  VL.  (Note  that  for  simplicity  we  are 
doing  our  analysis  in  the  tangent  plane  of  the  surface,  whereas  the  image  records  the  projection  of 
scene  radiance  from  the  surface  down  to  the  x-y  plane.  Thus,  our  analysis  ignores  the  effects  of 
projections.) 
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In  the  following  derivation,  the  operator  V  is  understood  to  be  the  gradient  taken  in  the  tangent 
plane,  or  u-v  plane,  at  a  point  on  the  surface.  A  short  calculation  shows  that  the  gradient  of  the 
scene  radiance 

VL  =  a(s  ■  V)h  +  (Va)s  ■  h.  (8) 

Eq.  8  highlights  two  factors  that  determine  the  gradient  of  the  scene  radiance.  The  first  term,  called 
the  geometric  gradient,  is  the  contribution  from  geometric  changes;  the  second  term,  called  the 
reflectance  gradient,  is  the  contribution  from  changes  in  the  BRDF. 

Consider  the  geometric  gradient  term: 

(s  ■  V)fl  =  UK,uSu  +  VK,vSv  (9) 


where  ku  and  kv  are  the  two  principal  curvatures,  and  su  and  sv  are  the  u  and  v  components  of  the 
light  source  in  the  u-v-n  coordinate  system. 

Consider  next  the  reflectance  gradient  term: 


->  da  „  da 
Va  =  u—  +v— 
du  dv 


(10) 


where,  after  a  little  calculation, 


da 

du 


da 

dv 


uc  _±  as  „  _=> 

Vca  •  t-  +  •  —  +  u  ■  Va 

du  du 

kuv  ■  (v^a  x  e  T  V§a  x  sj  +u  ■  Va, 
Kvii  ■  (v^a  x  c  +  Vs«  x  s')  +  ij  ■  Va 


where  V c  is  the  gradient  taken  with  respect  to  c  and  Vs-  is  the  gradient  taken  with  respect  to  s. 
We  combine  the  geometric  and  reflectance  gradients  to  get  an  overall  expression  for  VL : 


VL  =  (u  kusu  +  v  kvsv)  +  (ku  u  v  +  kv  v  u)  •  (Vc«  x  c  +  V -sa  x  s)  +  Va  s  ■  h  .  (11) 

s  __  ✓ 

v  V 

geometric  gradient  reflectance  gradient 

(Note  that  uv  and  vu  are  tensor  products.  The  associative  rule  applies  here.)  The  first  term  in  the 
square  bracket,  although  classified  as  part  of  the  reflectance  gradient,  is  induced  by  the  rotation  of 


14 


the  u-v-n  coordinate  system  along  the  lines  of  curvature.  If  the  BRDF  is  Lambertian,  the  above 
expression  becomes 

VL  =  (u  nnsu  +  v  kvsv)  +  (Va)s  •  h  .  (12) 

-  - - - '  ' - V - ' 

geometric  reflectance 

Now  that  we  have  the  expression  for  the  gradient  of  the  scene  radiance  in  terms  of  the  differential 
geometric  and  reflectance  properties  of  the  patch  of  surface,  we  can  determine  the  distribution  for 
scene  radiance  by  imposing  a  distribution  on  light  sources.  We  consider  light  sources  above  the 
tangent  plane  satisfying  sn  >  0.  We  felt  it  important  for  this  discussion  to  choose  a  distribution 
that  does  not  favor  any  particular  direction,  or  a  distribution  that  depends  on  the  magnitude  |s| 
alone,  even  though  it  may  prove  useful  to  assume  otherwise  under  some  particular  circumstances. 
In  addition,  we  assume  the  components  of  the  source  su.  sv,  and  sn  are  distributed  independently. 
With  these  assumptions,  it  can  be  shown  that  the  probability  density  for  s  is  given  by 


Ps(s)  =  J-  e  2 £a(*S+*3+*n),  SnG[0,oo). 


(\/2Tra 


(13) 


Note  again  this  distribution  is  uniform  over  the  direction  of  the  source.  We  subsequently  deduce 
the  corresponding  probability  density  function  of  VL  as  follows: 

General  Case: 

The  geometry  and  reflectance  factors  are  mixed.  The  distribution  in  u-v  coordinate  is 


Pu,v  (ll,  v) 

l  r°° 


Ps( 


JO 
1 

7T2  azKuKvC 


u  —  asr 


+ 


v  -  bsn  \  2  2 


+  sn)  d sn 


J_a£±K 

yplcrc 


e  T  dr. 


where  a  = 


V 

Ky 


(14) 


>  6  =  t  •  <=  =  \fW+W+~1’  t  =  Z- '  “d  c 

We  pay  particular  attention  to  two  special  cases  which  we  feel  are  the  determining  factors  for 


the  distribution  of  the  image  gradient. 


Special  Case  I:  Non-zero  Curvatures  with  Constant  Albedo 


15 


radius 


Figure  1:  Distribution  pu  ,v  for  constant  albedo  as  described  in  Eq.  15  in  polar  coordinates  (r,  ip). 


If  the  surface  patch  has  spatially  homogeneous  reflectance  (constant  albedo),  then  =  0.  In 

coordinate  system  u-v-n,  the  probability  density  function  for  VE  is  then 

7T  2  (J2KUKV 

Note  that  the  level  curves  of  this  function  are  concentric  ellipses  and  that  there  is  a  ridge  along  the 
larger  of  the  curvature  directions.  In  polar  coordinates  (r,  p),  where  r  =  y/u2  +  v2  and  tan  p  = 
there  are  two  equal  ridges,  at  p  =  0  and  ip  =  n,  along  the  direction  of  r-axis  whenever  44  7^  1 
as  shown  in  Figure  1.  The  ridges  grow  sharper  as  the  ratio  deviates  farther  away  from  1.  The 
existence  of  these  ridges  implies  that  the  gradient  is  more  likely  to  point  in  the  larger  principal 
curvature  direction.  As  stated  previously,  if  one  of  the  principal  curvatures,  say  k,v,  is  0,  such  as  in 
the  case  of  a  cylindrical  surface,  the  image  gradient  always  points  along  the  direction  of  the  nonzero 
curvature  k,.u. 


(15) 


Special  Case  II:  Zero  Curvatures  with  Nonconstant  Abedo 

Let  ku  =  k,v  =  0.  The  two  curvatures  are  the  same,  so  the  orientation  of  the  u-v  coordinate  system 
is  ambiguous.  We  may  simply  choose  the  direction  of  the  albedo  gradient  as  the  positive  u  axis. 
The  probability  density  function  for  the  gradient  is  then  simply 


Pu,v  {u,v) 


J2 _ L/JL^2 

—  e  S(v)0(u), 

707 


(16) 
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Figure  2:  Distribution  pu,v  for  a  planar  surface  as  described  in  Eq.  16  in  polar  coordinates  (r,ip). 

where  au  =  S  is  the  Dirac  delta  function,  and  ©  is  the  Heaviside  theta  function  (or  step  function). 
The  distribution  in  the  u-v  coordinate  system  is  one  thin  slice  along  the  positive  u  axis.  In  polar 
coordinates  (r,tp),  there  is  only  one  slice  at  <p  =  0  along  r  as  depicted  in  Figure  2.  Therefore,  the 
image  gradient  is  always  directed  along  the  albedo  gradient,  regardless  of  the  direction  of  the  light 
source. 

We  surmise  that  the  above  two  special  cases  dominate  the  natural  occurrence  of  images.  In 
most  natural  images,  the  gradient  in  intensity  will  be  largely  due  to  albedo  change  (case  II),  and 
the  geometric  influence  comes  into  play  only  when  the  material  is  almost  optically  homogeneous 
(case  I).  We  conjecture  that  the  linear  combination  of  these  two  cases  is  enough  to  describe  the 
distribution  of  the  image  gradient. 

We  want  to  develop  an  illumination  insensitive  measure  for  image  comparison  encompassing 
the  two  determining  factors  described  above.  If  we  are  given  the  directions  and  magnitudes  of 
principal  curvature  along  with  the  variance  of  the  BRDF,  then  we  can  use  the  density  function  pu,0 
to  determine,  in  a  probabilistic  sense,  how  faithful  any  image  is  to  the  given  values.  The  problem  we 
are  interested  in  is  slightly  different:  we  want  to  compare  two  images  and  determine  the  likelihood 
that  they  have  been  produced  by  the  same  object.  For  this  problem,  the  magnitudes  and  directions 
of  surface  curvature  and  the  BRDF  are  unknown.  We  must  then  look  at  the  joint  density  for 
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gradients  of  scene  radiance  as  given  in  two  images,  (as  they  are  the  only  observables  amongst  these 
previously  introduced  quantities)  and  integrate  out  the  unknown  unobservable  quantities,  i.e.,  the 
magnitude  and  direction  of  the  principal  curvatures  and  the  reflectance  parameters. 

Furthermore,  we  do  not  know  the  variation  in  the  strength  of  the  source  (the  standard  deviation 
a  in  Eq.  13).  This  variation  in  the  strength  of  the  light  source  is,  of  course,  embedded  in  —  but 
not  determining  —  the  variation  in  the  magnitude  of  the  gradient.  Because  of  this,  we  elect  to 
integrate  out  the  magnitude  of  the  image  gradient  as  well,  distilling  the  discriminative  power  of 
the  distribution.  On  the  other  hand,  the  magnitude  of  the  gradient  is  not  completely  proportional 
to  the  power  of  the  light  source.  As  the  image  intensity  decreases  as  in  the  shadow  region,  the 
variation  in  angle  ip  of  Pu,v{r,ip)  decreases,  as  can  be  deduced  from  Eq.  14,  and  is  clearly  depicted 
in  Figure  1.  Moreover,  the  signal  to  noise  ratio  decreases  as  well.  There  is  thus  distinction  between 
pu,v  at  high  image  intensity  and  low  image  intensity.  We  loose  some  information  in  collapsing  the 
dimension  of  the  probability  density  pv,v;  We  believe  nevertheless  that  under  usual  circumstances 
the  information  thus  lost  is  minute  and  we  gain  much  in  the  resulting  simplicity  and  computational 
speed.  The  experimental  result  in  Section  5  confirms  our  intuition. 

The  gradients  are  observed  in  a  fixed  x-y  coordinate  system.  Given  the  angle  7  between  x  and 
•u,  the  principal  curvatures  denoted  by  k’s,  and  the  albedo  gradient  Vo,  the  probability  density  of 
observing  a  gradient  with  magnitude  r  and  angle  ip  from  x  is  pr{r,  (p\y, «,  Vo)  =  r  P(u.V)(r  cos (<7  — 
7),rsin(</?  —  7)).  (The  scaling  of  p  by  r  comes  from  the  Jacobian  in  converting  from  Cartesian 
to  polar  coordinates.)  Noting  that  the  angular  dependence  is  only  on  the  difference  of  the  angles 
ip  —  7,  we  rewrite  the  density  as  pr{r,  ip  —  7|k,  Vo).  Now,  the  joint  probability  density  of  observing 
two  scene  radiance  gradients  (77 ,  <p\)  and  (r2,  ip 2 )  under  two  independent  and  identically  distributed 
light  sources  is 


p{ri,ipi,r2,<P2)  = 


J  Prin , 


ipi  -  7| K,(3)pr{r2,  P2  ~  (3)  dP(y,  k,/3), 


(17) 


where  (3  =  P( 7,  k,  (3)  is  the  probability  measure  on  the  unobservable  random  variables,  and 

the  integration  is  over  the  whole  sample  space.  Furthermore  if  azimuthal  symmetry  holds  for  P , 
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then  the  density  p  above  can  be  rewritten  as  a  function  of  three  variables 


p(n.i<P  =  <Pi  ~  <P2,r2)  =  f  f  pr{n,  {(pi  -  <p2)  ~  'y\K,/3)pr{r%i-'y\K,f3)  dydPjre,/?). 

J  J  7=— 7T 


(18) 


Azimuthal  symmetry  is  intrinsic  for  any  set  of  reasonably  random  image  samples:  the  unrestrained 
relative  rotation  of  the  objects  and  the  camera  along  the  optical  axis  of  the  lens  will  almost  surely 
render  the  azimuthal  angle  indistinguishable.  Eq.  18  implies  that  the  joint  density  depends  on  the 
magnitude  of  the  angle  between  two  image  gradients  or  the  absolute  value  of  <p  =  <p2  —  <p\.  If  is  thus 
an  even  function  with  respect  to  ip. 

In  keeping  with  our  earlier  assumption  that  the  density  function  pu.u  is  the  linear  combination  of 
the  special  cases  I  and  II,  we  evaluate  the  joint  density  p{r\.  ip,  r*2 1 /c,  (3)  as  such  a  linear  combination. 
Let  pi  and  p2  denote  the  distribution  in  polar  coordinate  in  the  special  case  I  and  II,  respectively. 
The  joint  density  for  special  case  I  is  then 


pi{ri,ip,r2\K) 

/7T 

pi{ri,'y\K)pi{r2,'y  +  ip\n)  dy 

-7T 

*(■■;+•■?)  (y^) 


k2i  k\ 


rf  +  r‘\  +  2r2r%  cos (2 ip)  I  ,  (19) 


(7TCT2 K,,UK,V)2  V  (2c)2 

where  Jq  is  the  O’th  hyperbolic  Bessel  function.  The  joint  density  for  special  case  II  is 


P2{n,(p,r2)  = 


/: 


P2{ri,j\(3)pi{r2,j  +  (p\f3)  dy 


7T2<72 


e  :,'di  (i(r-). 


(20) 


The  joint  density  is  then  the  linear  combination  of  Eq.s  19  and  20, 


p(ri,</?,r2|K,/3)  =  wipi  +  w2p2, 


(21) 


where  w  1  >  0  and  w2  >  0  are  the  weights  for  the  two  cases,  and  w  1  +  w2  =  1. 

Figure  3(a)  shows  a  graph  of  the  joint  probability  density  function  p(r\,  ip,  r2\n,  (3).  For  this 
graph  we  chose  a  =  1.7,  ku  =  4,  kv  =  1,  au  =  6,  w\  =  0.22,  and  w2  =  0.78.  We  chose  these 
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values  simply  to  illustrate  the  nature  of  the  shape  of  the  density.  Different  values  will  of  course  yield 
different  densities,  but  their  underlying  shapes  remain  qualitatively  the  same. 

Since  the  magnitude  r  of  an  image  gradient  is  directly  proportional  to  the  power  of  the  light 
source,  and  we  consider  only  the  local  effect,  the  magnitude  is  greatly  influenced  by  the  light 
source  power  and  thus  not  as  sensitive  to  the  underlying  geometry  and  photometry  as  the  angle  ip. 
Therefore,  if  we  desire  an  even  simpler  illumination  insensitive  measure,  we  can  integrate  p  over  r\ 
and  v2  and  obtain 

P<p{v)  =  J  j  P(ri,<P,r2)  dndr2.  (22) 

Clearly,  p v  inherits  from  p(r\,ip,r2)  the  property  of  having  a  unique  maximum  at  ip  =  0.  On  the 
other  hand,  the  original  measure  p(r\,  ip,  r2)  may  have  more  discriminatory  power,  especially  for  the 
dark  regions  of  the  images,  or  the  small  magnitudes  r\  and  r2. 

The  joint  distributions  pv(ip)  and  p(r\,ip,r2)  are  the  illumination  insensitive  measures  we  seek. 
Under  each  of  these  distributions,  we  would  expect  high  probability  to  be  assigned  to  two  images  of 
the  same  object  (differing  only  in  illumination,  not  viewpoint)  and  low  probability  assigned  to  two 
images  of  different  objects.  Yet,  we  are  hampered  by  our  ignorance  of  the  probability  distribution 
for  the  unobservables  (curvatures  and  albedo  gradients)  and  thus,  cannot  perform  the  integration 
in  Eq.  18.  We  note  that  all  of  the  above  analysis  is  done  in  the  tangent  plane  of  the  surface,  and 
consequently  ignores  the  effect  of  the  projection  of  the  radiance  from  the  surface  onto  the  image 
plane.  However,  it  will  become  apparent  in  the  next  section  that  this  effect  is  small. 

To  circumvent  this  difficulty,  and  to  offer  a  practical  solution,  in  the  next  section  we  will  construct 
empirically  the  distribution  in  Eq.  18  from  real  images  of  objects  under  varying  illumination.  We 
will  also  compare  the  empirical  data  against  the  theoretical  construct. 

5  Empirical  Joint  Density 

Using  a  geodesic  dome  (see  Figure  4)  with  64  photographic  flashes,  we  gathered  a  database  of  1,280 
images  of  20  objects  (see  Figure  5).  The  64  flashes  are  positioned  on  the  dome  to  cover  slightly  more 
than  a  hemisphere  of  directions.  The  objects  included  folded  cloth,  a  computer  keyboard,  cups,  an 
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Figure  3:  (a)  Joint  probability  densities  calculated  from  the  theoretical  model  of  the  two  image 
gradients  p(ri,ip,r2  =  50),  expressed  as  a  function  of  the  magnitude  of  one  gradient  and  the  angle 
between  the  two,  with  the  other’s  magnitude  set  to  50.  (b)  Empirical  data  of  the  same  function  as 
in  (a). 


umbrella,  plants,  a  styrofoam  mannequin,  among  others2.  A  stationary  camera,  positioned  at  the 
center  of  the  hemisphere  covered  by  the  flashed  light,  captured  64  images,  each  illuminated  by  one  of 
the  64  lights  flashing  in  quick  succession,  for  each  stationary  object.  The  64  images  of  a  telephone, 
one  of  the  objects,  is  shown  in  Figure  6.  We  estimate  the  density  p(ri,cf>,r2)  in  Eq.  18  directly 
from  a  histogram  of  the  image  gradients.  Eq.  18  is  used  with  the  except  that  the  integral  becomes 
discrete  summation  and  the  summation  is  taken  over  the  pixel  coordinates  and  objects  instead  of 
the  curvatures  and  albedo  gradient.  Since  the  image  function  is  scaled  by  the  power  of  the  light 
source,  we  only  need  to  collect  the  images  under  light  sources  with  the  same  power,  which  is  the 
case  for  the  flashs  on  the  dome.  The  image  gradient  distribution  under  light  sources  with  spherically 
symmetric  distribution  is  then  attained  by  integrating  the  distribution  p(r\,  </>,  7^2  |-Sl  =  1,  §2  =  1) 
under  unit  light  source  power  with  respect  to  the  power  distribution 

p{n,(p,r2)  =  [  f  p{  —  ,<p,  —  \si  =  1,52  =  l)e^{Sl+s'2)s\sl  dsid.s2,  (23) 

2707°  Jo  Jo  Si  S 2 

As  we  have  discussed  in  Section  4,  if  convenience  or  simplicity  is  desired,  p^ip1)  is  another 

2 All  databases  used  in  this  paper  are  available  for  download  from  subdirectories  “hrlfaces”,  “yaleAselected,”  and 
“objects”  at  ftp://Plucky.cs.yale.edu/FTPRoot. 
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Figure  4:  Geodesic  dome  with  64  flashes  used  to  capture  data  base  images. 


excellent  illumination  insensitive  measure.  The  distribution  p^(p)  in  angle  alone  can  be  gathered 
directly  from  the  image  database  under  a  uniform  single  power  lighting  distribution.  We  see  this  by 
substituting  the  above  expression  for  p(ri,ip,r2),  the  angular  distribution 


Pipip)  =  J  J  p(ri,ip,r2)  dr idr2 

/>(//!.  r-.  i)2h\  =  1,  s2  =  l)E[s]2dr]idr]2 
7ro^  ./  ./  e3°T(Sl+S2)<si's2  dsids2  j  J p{r]i,ip,r]2\si  =  1,  s2  =  l)dr]idr]2 


2n  a6 
2^2 

=  —  /  /  p{rji,  tp,  tj2\si  =  1,52  =  1)  drji dr/2 . 


(24) 


where  all  the  integrations  are  over  the  positive  real  axis,  and  E  is  the  expectation  value  with  respect 
to  the  radial  distribution  of  the  light  power  A.  We  perform  a  change  of  variables  from  (r,  A)  to  (r/.  A) 
to  obtain  the  second  equation.  As  shown  in  the  second  line,  Eq.  24  is  valid  for  any  radial  power 
distribution. 

A  slice  of  the  joint  probability  density  p(ri,(p,r2)  for  a  fixed  r2  is  shown  in  Figure  3(b).  As 
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Figure  5:  20  objects  in  the  database  used  to  collect  the  empirical  distribution  of  image  gradient. 


Figure  6:  A  telephone,  one  of  the  20  objects,  under  the  illumination  of  64  flashes. 


expected  in  Section  4,  the  distribution  exhibits  most  of  the  features  in  the  theoretical  model  depicted 
in  Figure  3(a).  It  is  symmetric  with  respect  to  the  plane  tp  =  0.  There  is  a  prominent  ridge  at 
p  =  0  and  it  is  the  unique  global  maximum  of  p(ri,(f>,r2)  on  the  line  of  arbitrarily  fixed  ri,r2- 
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As  anticipated  in  Section  4,  the  linear  combination  of  the  two  special  cases  captures  the  dominant 
properties  of  the  joint  image  gradient  distribution,  and  the  effect  of  projecting  the  tangent  plane 
to  the  imaging  plane  can  be  safely  ignored.  This  shows  that  the  statistical  regularity  of  the  scene 
radiance  gradient  does  reflect  the  intrinsic  geometric  and  reflectance  properties  of  surfaces  and  that 
this  regularity  can  be  exploited.  In  Section  6,  we  will  demonstrate  this  on  the  problem  of  face 
recognition  under  varying  illumination. 

6  Application  to  Face  Recognition 

Given  an  object  o  which  generates  images  I  and  J  under  two  different  lighting  conditions,  the  joint 
probability  of  observing  the  gradients  of  I  and  J  is  assumed  to  satisfy 

p(/,j)  =  IIMV4VJ,) 

■U  M 

=  n 

ieM 

where  M  is  the  set  of  pixel  indices  and  99  is  the  angle  between  the  two  gradient  vectors.  We  treat 
the  points  on  the  surface  as  being  independent,  ignoring  correlations  arising  from  spatial  proximity. 
From  a  Bayesian  perspective,  this  is  equivalent  to  assuming  that  when  the  images  come  from  two 
different  objects,  the  difference  in  the  gradient  direction  will  be  uniformly  distributed.  We  use 
density  function  p ^  as  the  illumination  insensitive  measure  for  the  sake  of  simplicity. 

We  apply  this  scheme  to  face  recognition.  Given  a  test  image  I  of  a  face,  we  compute  P(I,J)  for 
every  training  image  using  an  empirically  collected  probability  database  as  described  in  Section  5. 
The  one  training  image  having  the  highest  P  value  is  deemed  the  most  likely  to  have  come  from  the 
same  face  as  the  test  image  I.  Figure  7  shows  the  result  of  a  face  recognition  test  and  compares  it  to 
those  of  other  methods.  450  images  of  10  faces  each  under  45  different  lighting  conditions  are  used. 
One  image  of  each  face  under  frontal  illumination  is  taken  as  a  training  image.  The  recognition  test 
is  then  performed  for  the  remaining  440  images.  The  results  are  grouped  into  4  subsets  according  to 
the  lighting  angle  with  respect  to  the  frontal  or  camera  axis.  The  first  subset  covers  lighting  angles 
0°  —  12°,  the  second  covers  12°  —  25°,  the  third  25°  —  50°,  and  the  fourth  50°  —  77°.  We  compared 
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Comparison  of  Recognition  Methods 

Method 

Error  Rate  (%) 

vs.  Illumination 

Subset 

Subset 

Subset 

Subset 

1 

2 

3 

4 

Correlation 

0.0 

0.0 

11.7 

65.0 

Eigenfaces 

0.0 

0.0 

16.7 

69.3 

Eigenfaces 
w/o  1st  3 

0.0 

0.0 

3.3 

57.9 

Linear  subspace 

0.0 

0.0 

1.7 

12.9 

Cones-attached 

0.0 

0.0 

0.8 

9.3 

Gradient  Angle 

0.0 

0.0 

0.0 

1.4 

Cones-cast 

0.0 

0.0 

0.0 

0.0 

Figure  7:  Each  of  the  methods  except  Gradient  Angle,  which  has  only  one  training  image,  is  trained 
on  images  with  near  frontal  illumination  (Subsets  1).  This  graph  shows  the  error  rates  under  more 
extreme  light  source  conditions. 


our  method  with  those  tested  and  reported  in  [19].  The  image  gradient  method  clearly  outperforms 
all  other  methods  except  Cone  Cast. 

We  also  perform  the  test  using  the  probability  density  p(r\,  tp,  r2)  including  the  gradient  magni¬ 
tudes  instead  of  with  only  the  gradient  angle,  and  the  error  rates  are  on  a  par  with  those  above. 
It  confirms  our  earlier  assertion  in  Section  4  that  not  much  information  is  lost  when  the  magnitudes 
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Subset  4 


Figure  8:  Images  of  one  of  the  10  individuals  under  the  4  subsets  of  lighting. 

are  integrated  out.  Nevertheless,  we  hasten  to  reassert  that  when  the  image  intensity  is  extremely 
low  as  in  the  shadow  region  and  the  shadow  region  is  large,  we  expect  that  the  higher  dimensional 
density  would  be  more  discriminating. 

It  should  be  noted  that  all  the  methods  other  than  Gradient  Angle  use  all  of  Subset  1  for 
training,  and  so  by  definition  have  zero  error  rates  for  Subset  1.  In  contrast,  our  method  uses  only 
one  frontal  image  in  Subset  1  for  each  individual.  The  result  of  our  method  is  still  better  than  the 
others  with  the  exception  of  Cone  Cast.  Also  note  that  Cone  Cast  is  the  Illumination  Cone  method 
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in  [18]  in  which  Subset  1  images  are  used  to  construct  person  specific  illumination  representations. 
Our  method,  however,  has  the  advantage  of  being  much  simpler  and  faster  than  Cone  Cast  using- 
only  local  image  comparisons. 

It  is  worthwhile  to  note  that  the  probability  distribution  used  to  perform  the  test  is  gathered 
from  images  of  objects  rather  than  human  faces.  The  distributions  collected  from  different  categories 
of  objects  or  faces  are  remarkably  similar.  This  is  expected  from  the  analysis  in  Section  4. 

7  Conclusion 

This  paper  presents  two  results:  first,  illumination  invariants  do  not  exist  for  Lambertian  surfaces; 
and  second,  the  angle  (or  direction)  of  the  image  gradient  is  insensitive  to  changes  in  illumination 
direction.  The  latter  statement  is  consistent  with  the  conclusion  in  [1]  that  linear  filters  for  image 
comparison  do  not  exist,  since  the  gradient  angle  is  a  nonlinear  function  of  the  image.  However,  we 
cannot  conclude  that  image  edges  are  good  measures  of  image  comparison  under  varying  illumination 
—  in  fact  the  contrary  is  true.  Most  edge  detection  methods  are  highly  sensitive  to  the  magnitude 
of  the  image  gradient.  As  we  can  see  from  Eq.  12,  the  magnitude  of  the  gradient  varies  drastically 
with  the  change  in  the  direction  of  the  light  source.  The  distributions  in  Figure  3  also  shows  the 
slow  variation  of  the  density  function  with  respect  to  the  magnitude  of  the  gradient.  Therefore 
the  gradient  magnitude  is  a  poor  indicator  of  the  underlying  surface  geometry  and  photometry. 
Nevertheless,  when  combined  with  the  gradient  angle,  the  magnitude  may  extract  more  information 
from  extremely  low  intensity  images. 

Our  assumption  about  pixel  independence  implicit  in  Eq.  25  is  obvious  incorrect  —  the  light 
sources  are  fixed  for  all  pixels  in  a  given  image,  and  neighboring  surface  points  tend  to  have  sim¬ 
ilar  geometric  and  reflectance  properties.  We  hope  to  remedy  the  crudeness  of  our  independence 
assumption  by  exploring  how  spatial  dependencies  in  pixel  values  [2]  can  be  exploited  to  produce 
even  more  discriminating  measures  of  the  images. 
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8  Appendix 

The  propositions  in  Section  3  are  proved  here. 

The  following  observation  on  the  characteristic  curve  Eq.  5  is  crucial  to  the  proofs  of  the  lemmas. 

Observation  8.1  Let  q  be  an  arbitrary  point  in  i?3,  Sq  the  plane  containing  q  parallel  to  l  and  s,  and 
rq  the  characteristic  curve  passing  through  q.  It  is  clear  that  if  fq  exists,  it  remains  in  the  plane  Sq. 
Moreover,  the  characteristic  vector  at  every  point  resides  in  the  planar  cone  C  =  {al  —  bs,  a,  b  >  0} 
generated  by  l,  —s,  and  the  negative  of  C.  Consequently  the  curve  rq  is  also  confined  in  the  two-sided 
cone  q  +  C  U  q  —  C .  Specifically,  the  projection  pp  of  rq  onto  the  XY-plane  is  contained  in  a  cone 
Cp,  the  projection  of  C  on  the  XY-plane. 

We  also  need  another  lemma,  Whitney’s  Extension  Theorem  [15],  to  extend  a  C1  function  onto 
a  larger  set  in  the  proof  of  the  lemmas  in  the  main  text. 

Lemma  8.1  Let  A  be  a  closed  set  of  points  a  in  Rm  at  which  the  values  and  derivatives  of  a  desired 
C1  function  are  prescribed  by  linear  polynomials  Pa  :  Rm  — >  R.  For  each  compact  subset  C  of  A 
and  6  >  0,  let  f(C,S)  be  the  suprem.um.  of  the  numbers  \Pa  '^1^'  ,  ||  DPa(b)  —  DPi{b)\\,  where  D 

is  the  differentiation  operator,  over  all  a,b  <E  C  with  0  <  |a  —  b\  <  S.  If  the  prescribed  data  satisfy 
the  coherence  condition  that  lim^o  £(C,  S)  =  0  for  each  compact  subset  C  of  A,  then  there  exists  a 
C 1  funciton  g  satisfying  g{a )  =  Pa{a),  Dg{a)  =  DPa{a ),  Va  G  A,  and  mfxeR™  g{x)  =  inf aeAPa{a), 
SUPxe-Y  =  suPae.4Pa(«)- 

Proof  of  Lemma  3.3. 

We  shall  first  caution  the  reader  that  this  lemma  concerns  the  global  existence  of  a  solution  as 
opposed  to  the  local  existence.  We  intend  to  turn  the  ordinary  equation  into  an  integral  equation 
and  solve  the  integral  equation  with  iterations.  However,  the  interval  of  integration  should  be  fixed 
and  the  trajectory  from  each  iteration  should  all  reach  the  boundary  of  the  domain  f l.  Thus  the 
vector  field  should  be  extended  to  a  larger  domain  containing  f 1  to  satisfy  these  requirements. 
Proof.  Without  loss  of  generality,  we  assume  that  the  given  point  is  at  the  origin  and  f l  is  bounded 
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by  R  >  0.  Since  I  and  J  are  positive  and  continuous  on  a  compact  set,  M'  >  /,  J  >  c  for  some 
constants  M'  >  c  >  0.  Recall  (Observation  8.1)  that  the  characteristic  vectors  in  the  XY-plane  are 
confined  in  the  cone  Cp.  We  denote  by  h  the  vector  bisecting  the  angle  subtended  by  lp  and  —  sp.  As 
I  and  J  are  both  bounded  from  below  by  c  >  0,  the  projection  of  the  vector  field  on  h  is  bounded 
from  below  by  c(lp  +  sp)  cos  6,  where  6  is  half  the  angle  between  lp  and  sp,  the  projections  of  l  and 
s  on  the  XY  plane  respectively. 

If  9  <  j,  any  curve  starting  from  the  origin  with  these  velocity  (or  tangential)  vectors  will  leave 
Q  in  time  T  =  cc^sff-  If  0  =  j,  Ilp  —  Jsp  is  on  a  straight  line.  The  exit  time  is  then  simply  bounded 
by  T  =  Let  max^Q  \I(p)lp  —  J(p)sp\  =  M.  By  setting  Pfi(p)  =  I(a)  +  DI(a)(p),  and  the  same  for 
J(p),  Lemma  8.1  (Whitney’s  Extension  Theorem)  allows  us  to  extend  the  functions  I,  J  G  C1)!!)  to 
C1  functions,  again  denoted  by  I  and  J,  onto  the  closed  ball  B  of  radius  MT ,  such  that  both  newly 
constructed  I  and  J  maintain  their  original  values  on  il.  and  assume  the  same  respective  maxima 
and  minima  on  B  as  on  Q. 

We  are  now  in  a  position  to  recursively  construct  a  sequence  of  functions  {pn}^L0  :  [0,  T]  — >  B. 
Let 

Po{t)  =  0.  (25) 

Suppose  pn  is  defined  on  [0,T],  Define 

Pn+i{t)  =  [l{pn{T))lp  -  J(pn(r))sp]  dr,  t  6  [0,  T].  (26) 

Clearly  pn(t)  E  B,  Vt  E  [0,  T],  Vn  E  N  or  {p}%L o  is  uniformly  bounded  by  B.  Since  I,  J  E  C1 
on  the  convex  and  compact  set  B ,  by  the  mean  value  theorem,  \(I(p2)lp  —  J(p‘>)f?p)  —  (I(p\)lp  — 
J(pi)sp)\  <  \\p2  —  pj|,  V/l| ,  p‘>  E  B  for  some  constant  A  >  0.  This  leads  to 

\pn+i(t)  ~  Pn(t) I  <  A  /  I Pn(r)  -  pn- i(t)|  dr.  (27) 

Jo 

Apply  the  bound  \Il  —  Js|  <  M  to  Inequality  27, 

\pi(t)\  <  Mt.  (28) 
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We  obtain  by  induction, 


^n^n+l 

\pn+i(t)  -pn(t) I  <  M  ,  Vi  G  [0,  i]. 

(n  +  1)! 


(29) 


The  series  ^(pn+i(i)  —  Pn(t ))  is  majorized  by  that  of  jeA(.  Therefore  the  sequence  pn  tends 

n— 0 

uniformly  to  a  continuous  function  p.  Since  is  closed,  p(t)  G  _B,  Vi  G  [0,T]. 

lip  —  is  continuous  on  the  compact  set  B ,  and  is  therefore  uniformly  continuous.  Together 
with  the  uniform  convergence  of  {pn}n =o,  it-  enables  us  to  take  the  limit  n  — >  oo  on  both  sides  of 
Eq.  26,  and  substitute  pn  with  p  in  the  integrand. 

Let  (  =  infte[0T]{i  :  p(t)  G  B  —  Int( fl)}.  Since  the  image  of  the  continuous  function  p  on  the 
compact  set  [0,  T]  is  also  compact  and  therefore  closed,  its  intersection  with  the  closed  set  B  —  Int(£l) 
is  closed.  Again  by  continuity  of  p  {t  :  p(t)  G  B  —  Int(il)}  is  closed,  and  thus  p{()  G  B  —  Int(Q). 
Consequently  p(Q  G  8Q  and  p(t)  G  Int{ fi),  Vt  G  [0,  £).  Thus  p(t ),  t  G  [0,  (j]  is  a  solution  of  Eq.  5 
on  Q. 

The  uniqueness  of  the  characteristic  curve  through  the  origin  on  can  be  shown  following  the 
classical  argument.  Suppose  there  are  two  solutions  to  Eq.  5,  p\  and  f>2  with  the  same  initial  points 
pj(0)  =  p2 (0)  =  0  (They  may  exist  only  in  the  interior  of  without  reaching  the  boundary.). 
Suppose  they  are  both  defined  up  to  t.  We  have 

\p2{t)  -  pi{t)\  <  I  (I[p2(t)]Ip  —  J[P2{t)\s'p)  —  (I[pi(t)]Ip  —  J[pi(r)]5p)  dr 
Jo 

<  A  /  |p2(r)  -pi(r)|  dr, 

Jo 

or,  what  is  equivalent, 

^  lo  dT)  -°‘  (30) 

It  is  obvious  from  the  initial  condition  | p2 ( 0 )  —  /5j(0)|  =  0,  that  \p2{t)  —  pi{t)\  =0  for  every  defined  t. 
Consequently,  there  is  a  unique  solution  emanating  from  the  origin  and  evolving  positively  in  time. 
It  lands  on  the  boundary  at  t  =  (■ 

As  for  the  other  part  of  the  curve  in  the  negative  direction,  replace  t  with  —t,  every  step  of  the 
above  argument  carries  through.  □ 
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Proof  of  Lemma  3.4.  The  proof  for  these  global  properties  is  the  same  as  that  for  the  corre¬ 
sponding  local  ones  in  many  classic  texts  (e.g.,  [11]).  □ 

Proof  of  Lemma  3.5. 

Lemma  3.3  guarantees  that  every  point  in  0  is  an  outgrowth  from  the  boundary.  A  surface  covering 
the  whole  domain  O  can  be  constructed  by  growing  the  characteristic  curves  from  boundary  curves. 
However  the  boundary  conditions  should  not  lead  to  multi-valued  function  /.  To  avoid  this  difficulty, 
we  simply  expend  the  vector  field  to  a  triangle  containing  Q,  and  assign  boundary  values  to  one  side 
of  the  triangle  which  will  not  lead  to  multi-valuedness.  The  required  surface  is  simply  the  restriction 
of  the  surface  constructed  on  the  triangle  to  0. 

Proof.  Let  c  =  lp  —  sp).  Draw  a  straight  line  T  perpendicular  to  c  such  that  is  on  the  side  of 
T  that  c  points  to.  Draw  two  lines  parallel  to  lp  and  —  sp  respectively,  so  that  together  with  T  they 
form  a  triangle  enclosing  the  compact  set  0.  Then  restrict  T  to  denote  the  closed  segment  that  is 
a  side  of  the  triangle  (If  lp  and  —  sp  are  parallel,  we  can  either  consider  the  triangle  as  having  one 
apex  at  infinity,  or  cut  off  the  two  parallel  sides  and  form  a  rectangle,  so  long  as  it  contains  lb).  We 
again  use  Whitney’s  Extension  Theorem,  Lemma  8.1,  and  extend  I  and  J  on  D  to  two  C]  functions 
on  the  closed  triangle,  retaining  their  original  minima  and  maxima. 

Let  the  compact  set  D  in  the  premise  of  Lemma,  3.3  be  the  triangle  constructed  above.  It  then 
follows  from  Lemma  3.3  and  Observation  8.1,  and  the  fact  that  T  is  not  parallel  to  any  characteristic 
direction  of  the  PDE,  that  through  any  point  p  E  D  runs  a  unique  characteristic  line  emanating 
from  a  unique  point  on  T  and  reaching  the  other  two  sides  of  the  triangle.  Parameterize  T  by  its 
length  s  from  one  of  its  two  ends.  Denote  the  characteristic  curve  emanating  from  T  by  p(s,t).  By 
Lemma  3.4,  (x(s,t),y(s,t))  =  p{s,t)  establishes  a  C1  diffeomorphic  coordinate  chart  for  the  triangle 
(a  two  dimensional  manifold). 
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Let  fo{s)  be  an  arbitrary  C 1  function  on  T  and  let 

f(s,t)=  f  (I[p(s,  t )}lz  —  J[p(s,  t )]s2)  dr  (31) 

Jo 

It  is  obvious  that  f{s,t)  E  C1,  and  therefore  r(s,t)  =  (x(s,  t),  y(s,  t),  f(s,  t))  is  a  C1  manifold.  By 
the  diffeomorphism  (x,y)  =  p(s,t ),  /  could  be  written  as  a  function  of  (x,y)  and  the  manifold  is 
the  graph  of  /.  The  desired  surface  is  the  manifold  restricted  to  f l.  □ 


Proof  of  Lemma  3.2. 

As  is  mentioned  in  Section  3,  for  Eq.’s  2  and  3  to  adequately  provide  the  global  description  of 
the  required  surface,  it  is  necessary  that  the  surface  constructed  should  have  no  attached  or  cast 
shadows.  Otherwise  Eq.  1  instead  should  be  used.  Since  all  the  characteristic  curves  stay  confined 
in  a  planar  cone  described  in  Observation  8.1,  no  point  on  the  surface  will  attach  a  shadow  to  itself 
or  cast  one  on  another  point  so  long  as  the  boundary  curve  assigned  to  the  side  of  the  triangle 
satisfies  certain  condition. 

Proof.  Construct  the  same  triangle  as  in  Lemma,  3.5.  Let  fo{s)  be  the  initial  condition  on  T  as 
defined  in  the  proof  of  Lemma  3.5.  By  Lemma  3.5,  the  graph  can  also  be  denoted  as  f(s,  t  =  0).  In 
addition,  r(s,  0)  is  such  that 


>o. 

os 


(32) 


Point  a  casts  a  shadow  on  point  b  if  and  only  if  a  and  b  are  on  the  same  line  parallel  to  either  s  or 
l.  By  the  way  we  construct  the  surface,  it  is  a  manifold  parameterized  by  (s,  t)  via  a  diffeomorphism 
r(s,t).  Let  a  =  r(si,tfi)  and  b  =  r(s2,t2)1  we  distinguish  two  cases:  1)  .S|  =  ts2;  2)  .sj  /  s2. 

In  case  of  1), 


r(si,t2)  -  r(si,ti)  =  £  (l[p(si,t)]l  -  J[p(si,i)]s)  dr. 


(33) 


Since  J[p(s,t)]  >  0  and  continuous,  —  fj*  J[p{si,t)]  dr  <  0  and  thus  r(si,t2)  —  r{si,ti)  is  not 
parallel  to  l.  By  the  same  token,  r(.si,t2)  —  r(.si,ti)  is  not  parallel  to  s  either. 
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For  case  2),  we  show  that  r(si,  t2)  —  r(si,  t\ )  is  not  in  the  plane  of  l  and  s. 


-sxl-r(s,t)  =  — s  x  l 


uo 


(l[p(5,  t)]1  -  J[p(s ,  r)]s)  dr  +  f(s,  0) 
=  —s  x  l  ■  ^a(s,  i)i  —  b(s,  t)s  +  r(s,  0)  j 
=  — s  x  l  ■  f(s ,  0). 


(34) 


a(s,t)  and  b(s,t)  stand  for  the  scalar  functions  of  vectors  l  and  s.  Without  loss  of  generality,  we 
assume  s 2  >  S] . 


s  x  l  (f(s2,t2)  -  r(si,ti))  =  -sxl  ■  (f(s 2,0)  -  f(si,0))  >  0, 


(35) 


where  the  last  inequality  comes  from  integrating  Inequality  32  over  s  from  si  to  52  and  the  claim  is 
verified. 

Divide  Inequality  34  by  s2  —  si  and  let  s2  — >  -S| ,  we  obtain 


_s-xr.®^i>o. 

os 

Then  for  the  surface  normal  vector  x  dl  ,  the  sign  of  the  inner  product  is  always 


(36) 


dr(s,t)  dr(s,t)  _  dr(s,t) 
■  s  = 


ds 


dt 


ds 


(. it -  Js) 


x  s  >  0, 


(37) 


since  I(x,  y)  and  J(x,  y)  are  positive.  The  same  inequality  holds  for  the  inner  product  of  the  normal 
vector  and  l.  In  other  words,  there  is  no  attached  shadow.  □ 


Proof  of  Theorem  3.2. 

We  have  shown  that  in  the  absence  of  interreflections  we  can  construct  an  object  and  two  light 
sources  that  produce  the  two  images  I  and  J.  We  now  show  that  with  interreflections,  there  will 
still  be  no  continuous  discriminative  invariants.  To  do  this  we  show  that  if  we  scale  up  the  magnitude 
of  the  light  sources,  while  scaling  down  the  albedos  by  an  inverse  amount,  we  can  suppress  the  effect 
of  interreflections  to  the  infinitesimal  and  thus  approximate  the  given  images  with  those  generated 
by  the  surface  with  arbitrary  precision. 
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Proof.  An  image  I{x,y)  generated  by  an  interreflecting  surface  f{x,y)  on  a  compact  set  (x,y)  E 
can  be  described  by  the  differo-integral  equation 

I  =  Jo  +  aS[f]I,  (38) 


where  I  is  the  image  function,  and  Iq  is  the  intensity  of  the  first  reflection  off  the  surface  of  the 
light  directly  from  the  light  source.  Linear  operator  S[f]  is  a  function  of  the  continuous  surface 
/  E  C(fl)  fl  piecewise  C1^), 

smp)  =  -  [  d y  (39) 

Jn  rz 

where  p  =  (x,y),  p  '  =  (x',y'),  p,p'  E  ft.  d2p'  stands  for  area  element  on  the  surface,  f  = 
(x1  ~  x,y'  -  y,z'  -  z),  z  =  f(p)  =  f(x,y),  and  z'  =  f(p')  =  f(x',y').  r  =  \r\,  r  =  r/r.  h  and 
h'  are  the  surface  normal  at  points  (x,  y,  f(x,  y))  and  (a/,  y1 ,  f(x',  y')),  respectively.  5(p,p')  is  an 
index  function  equal  to  1  if  the  line  segment  connecting  p  and  p  '  does  not  intersect  the  surface 
(x,y,  f(x,y)),  and  0  if  the  line  segment  does. 

It  is  clear  by  conservation  of  energy,  the  maximum  norm  of  the  operator  ||S'[/]||l°°  <  1,  and 
0  <  arn  =f  max^gQ  a(p)  <  1.  The  solution  of  Eq.  38  is  then  readily  obtained  by  Neumann  series, 
i.e., 

OO 

/=£  anSnI0,  (40) 

n=0 

where  Sn  denotes  the  n’th  operation  of  the  operator  S  on  Iq  and  «n  is  the  n’th  multiplication.  Thus 


1 1 -I, 


0  \\L° 


<  Otr, 


'0||L° 


1  0?7l||*S'||ioo 


(41) 


Given  two  image  functions  Iq  and  Jq  and  two  arbitrary  light  source  direction  s  and  l.  there 
are,  according  to  Theorem  3.1,  a  surface  /  and  an  albedo  a,  such  that  Iq  =  as  ■  n  =  a  as  ■  h  and 
Jq  =  al  ■  n  =  bal  ■  h  (and  there  is  no  cast  or  attached  shadow  on  /),  where  a  =f  arn \ s\ .  b  =f  arn | / 1 , 


and  a  =f 

Ctm 


Define  /  and  J,  as  at  the  beginning  of  the  proof,  to  be  the  images  generated  by  the  surface  / 
with  interreflection  under  s  and  l,  respectively.  Let  s  =f  |s|  =  ~r—  and  l  =f  |/|  =  I{am)  and 
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J(am)  are  then  functions  of  0  <  am  <  1.  Iq  and  Jo,  on  the  other  hand,  are  independent  of  am.  It 


then  follows  from  Inequality  41  that  I(am)  — >  Iq  and  J{am )  — >  Jo,  as  a,m  — >  0+. 

Suppose  /t  is  a  continuous  illumination  invariant.  By  the  definition  of  illumination  invariant, 
n(I(am))  =  n(J(am)),  since  I  and  J  are  generated  from  the  same  surface  and  albedo.  By  the 
continuity  of  /i,  /.t(/o)  =  limQ,m_>0+  u{I{am))  =  lim am_>o+ M  J(am))  =  M^o)-  Hence  the  invariant  /i 
is  nondiscriminative.  □ 
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